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ABSTRACT 



Context. MHD turbulence is known to exist in shearing boxes with either zero or nonzero net magnetic flux. However, the way 
turbulence survives in the zero-net-flux case is not explained by linear theory and appears as a purely numerical result that is not well 
understood. This type of turbulence is also related to the possibility of having a dynamo action in accretion discs, which may help to 
generate the large-scale magnetic field required by ejection processes. 

Aims. We look for a nonlinear mechanism able to explain the persistence of MHD turbulence in shearing boxes with zero net magnetic 
flux, and potentially leading to large-scale dynamo action. 

Methods. Spectral nonlinear simulations of the magnetorotational instability are shown to exhibit a large-scale axisymmetric magnetic 
field, maintained for a few orbits. The generation process of this field is investigated using the results of the simulations and an 
inhomogeneous linear approach. We show that quasilinear nonaxisymmetric waves may provide a positive back-reaction on the large- 
scale field when a weak inhomogeneous azimuthal field is present, explaining the behaviour of the simulations. We finally reproduce 
the dynamo cycles using a simple closure model summarising our linear results. 

Results. The mechanism by which turbulence is sustained in zero-net-flux shearing boxes is shown to be related to the existence of a 
large-scale azimuthal field, surviving for several orbits. In particular, it is shown that MHD turbulence in shearing boxes can be seen 
as a dynamo process coupled to a magnetorotational-type instability. 
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1. Introduction 

The problem of angular momentu m transport is a central issu e 
of accretion disc theory. Following Shakur a~& Sunvaevl (Il973l) . 
angular momentum transport is often modelled assuming the 
disc is turbulent, using a kind of turbulent viscosity (the so- 
called a model). However, the way discs may become turbulent 
is still a highly debated subject. A possible r oute to turbulence 
is the magnetorotational instability (MRI, Balbus ^fe Hawlevl 
1991) which appears in discs sufficiently ionis ed to be coupled 
with the magnetic field lines dGammiej|1996l) . This instability 
has been extensively studied in its nonlin ear regime using lo- 
cal dHawlev et al.lll995t IStone et al.|[l996h and global dHawlevI 
2000) numerical simulations. More recently, the effect of non- 
ideal MHD has been investigated n umerically in shearing boxes 
with either zero {Fromang et alj2007h or nonzero ( Flemin g et alj 
2000; iLesur & Longarettill2007l) net magnetic flux, showing a 
strong dependence on the magnetic Prandtl number Pm. 

These new results bring back the question of the efficiency of 
MHD turbulence in accretion discs since Pm is believed to vary 
by se veral orders of magnitude in real objects dBalbus & Henril 
2008). In particular, as turbulence seems to disappear in zero- 
net-flux simulations for Pm < 1, the existence of a turbulent 
flow in low-Pm objects such as pr otoplanetary dis c s is qu estion- 
able. However, as pointed out by [Fromang et al.l d2007l) . today 
simulations can reach Reynolds numbers that are very small (up 
to 10 4 ) compared to those of real discs (~ 10 10 ). Therefore, no 
clear conclusion for accretion discs can be drawn from these nu- 
merical results. 

It has been point e d ou t by iPessah et"ai1 d2007l) and 
Fromang & Papaloizou ( 2007) that numerical resolution can 



also play an important role in simulations. In particular, increas- 
ing resolution in zero-net-flux simulations without any physi- 
cal dissipation leads to a weaker turbulence and transport. One 
might conclude from these results tha t turbulence should dis- 
appear in real astrophysical systems (Pes sah et al.l 12007). We 
think however that this conclusion is questionable since the ideal 
MHD model does not hold in turbulent flows, mainly because a 
dissipation scale necessarily exists, at which non-ideal effects are 
not negligible. Of course, ideal MHD simulations must include 
some sort of numerical dissipation in their algorithm, which 
then implicitly defines a dissipation scale of the order of the 
numerical grid scale. However, this algorithm-dependent dis- 
sipation is quite differen t from physical dissipation processes 
dLesaffre & Bal bus 2007), leading potentially to numerical ar- 
tifacts for large-scale properties, such as turbulent transport. 

On the other hand, the way MHD turbulence is sustained in 
these simulations is not well understood. It is known that when 
a mean vertical field is applied, the magnetorotational instability 
destabilises the flow and leads to developed three-dimen sional 
turbulence dBalbus & Hawlevl 1 1 99 U lHawlev et al]|l995l) . This 
picture is not directly applicable to zero-net-flux simulations 
since the magnetic field can be dissipa ted either by a finite re- 
sistivity, or by the turbulence itself (see lHawlev et al.l (Il996b for 
an extensive discussion of this point). Therefore, even if one as- 
sumes that the MRI appears locally because of a given magnetic 
field configuration, one has to regenerate this field with some 
sort of dynamo mechanism. The whole process sustaining turbu- 
lence in zero-net-flux simulations is therefore a nonlinear mech- 
anism, potentially involving a kind of magn etorotational insta- 
bility at some stage. As already mentioned bv lBalbus & Hawlevl 
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(1992), this means that this disc dynamo is intrinsically non- 
linear, as it requires a non-negligible Lorentz force, and can- 
not be described as a kinematic dynamo. Following these ideas, 
iRincon et al.l (|2007) obtained a steady nonlinear solution to the 
MHD equations in the zero-net-flux case, using rigid conducting 
walls as radial boundary conditions. This solution is clearly a 
first step toward the understanding of the mechanism working in 
shearing boxes, although the boundary conditions and Reynolds 
numbers are significa ntly different compared to turbulent simu- 
lations. Note also that lBranden burg et al] d 19951) studied this dy- 
namo process in discs, using boundary conditions allowing for 
mean flux variations. Although a azimuthal field was generated 
in their simulations, no physical understanding of the underlying 
process was provided. 

In this paper, we describe a possible mechanism able to sus- 
tain MHD turbulence in the shearing box with zero net flux. 
First, we recall the resistive MHD equations in the shearing box, 
and the numerical method used to solve them. Next, we inves- 
tigate the temporal evolution of some zero-net-flux simulations, 
and we exhibit a long-timescale cycle for the large-scale mag- 
netic field. The origin of this cycle is described using a spectral 
decomposition, and we demonstrate that its origin is related to 
some properties of the nonaxisymmetric structures of the flow. 
We then study a linear theory of nonaxisymmetric waves in 
the presence of a large-scale magnetic field similar to the one 
observed in simulations. We show that this linear theory pre- 
dicts the right properties for the nonaxisymmetric structures and 
partially explains the nonlinear cycle. We then provide a phe- 
nomenological closure model for the evolution of the large-scale 
field, based on our previous findings. We show that this model re- 
produces the basic behaviour of the cycles, and provides a mech- 
anism able to sustain turbulence in dissipative MHD shearing 
boxes. The existence of a long-timescale cycle and large-scale 
structures in these simulations is the most significant finding of 
our investigation, mainly because it shows that MHD turbulence 
is able to generate large scales, independent of the dissipative 
scales. This mechanism can also be seen as a way to generate 
large-scale magnetic fields, providing a dynamo mechanism spe- 
cific to accretion disc turbulence. These findings are discussed in 
detail in the concluding section, along with the implications for 
real astrophysical systems. 

2. Shearing-box equations and numerical method 

MRI-related turbulence has been extensively studied in the lit- 
erature. Therefore, we will recall here briefly the basic equa- 
tions for the sheari n g-box model. The r eader may consult 
lHawlev et ail (1 19951) . iBalbusI d2003l) and iRegev & Umurhanl 
(2008) for an extensive discussion of the properties and limi- 
tations of this model. Since MHD turbulence in discs is sub- 
sonic, we will work in the incompressible approximation, which 
allows us to eliminate sound waves or density waves. We also 
neglect vertic al stratification, consisten t with the local shearing- 
box model dRegev & Umurhanl 120 081 We include in our de- 
scription a molecular viscosity and resistivity to minimize the 
artifacts of numerical dissipation. 

The shearing-box equations are found by considering a 
Cartesian box centred at r — Rq, rotating with the disc at an- 
gular velocity Q = Q(/?o) and having dimensions (L x ,L y ,L z ) 
with Li «: Rq. We define Rq(/> — > x and r - Rq — * —y for con- 
sistency with the sta ndard notation for plane Couette flow (e.g. 
iDrazin & ReidllT98Tb . Note that this defini tion differs from the 
standard notation used in shearing boxes (lHawlev et al.| [l995) 
with x — > -ysB, y — > *sb and z — > Zsb- In this rotating frame, 



one eventually obtains the following set of equations: 



d,u + V • (u <g> u) = - Vn + V • (B ® B) 

-2fl X u + 2D.Sye y + vAu, (1) 

d,B = V X (u X JB) + t]AB, (2) 

V ■ u = 0, (3) 

V • B = 0. (4) 



The boundary conditions associated with this system are peri- 
odic in the x and z direction and shearing-periodic in the y di- 
rection (see lHawlev et al.l (1 19951) for a complete description of 
these boundary conditions). In these equations, we have defined 
the mean shear S = -rd r Q., which is set to S — (3/2)fi, assum- 
ing a Keplerian rotation profile. The generalized pressure term 
n includes both the kinematic pressure P/po and the magnetic 
pressure B 2 /2. One should note that the generalized pressure n 
is actually a Lagrange multiplier enforcing equation (O, and is 
therefore computed solving a Poisson equation. Note also that 
the magnetic field is expressed in Alfven-speed units, for sim- 
plicity. 

The steady-state solution to these equations is the local 
Keplerian profile u = Sye x . In this paper, we will consider only 
the turbulent deviations from this Keplerian profile. These may 
be written as v = u - Sye x , leading to the following equations 
for v. 

d,v + V • (v ® v) = - Vn + V • (B <g> B) - Syd x v 

+(2Q - S )v y e x - 2Qv x e y + vA v, (5) 
d,B = -Syd x B + SB y e x 

+ V x (v x B) + r]AB, (6) 

V • v = 0, (7) 

V • B = 0. (8) 

Following lHawlev et al.l d 1995b . one can integrate the induction 
equation (|6]l over the volume of the box, leading to: 

^-=S(B y ) ea> , (9) 
ot 

where () denotes a volume average. Therefore, the mean mag- 
netic field is conserved, provided that no mean radial field is 
present. This allows us to define the zero-net-flux shearing box, 
as the box in which (B) - 0. 

To numerically solve the shearing-box equations, we use 
a spectral Galerkin representation of equat ions I©-© in the 
sheared frame (see ILesur & Long aretti 2005). This frame al- 
lows us to use a Fourier decomposition since the shearing- 
sheet boundary conditions are transformed into perfectly pe- 
riodic boundary conditions. Moreover, this decomposition al- 
lows us to conserve magnetic flux to machine precision with- 
out any modification, which is an advantage compared to finite- 
difference or finite- volume methods (the total magnetic flux cre- 
ated during one simulation is typically 10~ n ). Equations © and 
© a re enforced to machine precision using a spectral projec- 
tion (Pevret 2002). The nonlinear terms are computed with a 
pseudospectral method, and aliasing is prevented using the 3/2 
rule. To always compute the physically relevant scales in the 
sheared fr ame, we use a remap met hod similar to the one de- 
scribed by Umurhan & Regevj d2004l) . This routine redefines the 
sheared frame every T Tem!lp = L X /(L V S) and we have checked 
that none of the behaviour we describe in this paper was related 
to this numerical timescale. Since spectral methods are very lit- 
tle dissipative by nature, we check that numerical dissipation is 
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Fig. 1. Time-evolution of B x , averaged in the x and y directions. Notice the long-lived (T ~ 50 S ') vertical structures, mostly found 
with long vertical wavelengths. 
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Fig. 2. Fourier analysis of B x (k x = 0, k y =0,k z = 2n/L z ): ampli- 
tude (top panel) and phase (middle panel). This mode exhibits 
long-timescale (T ~ 50 S -1 ) cycles during which the phase 
is approximately constant. The transport coefficient (a, bottom 
panel) does not show this cyclic behaviour. 



kept to very small v alues, computing the total energy budget at 
each time step (see iLesur & Long aretti 2005, for a discussion 
of this procedure). We therefore ensure that numerical dissipa- 
tion is responsible for less than 3% of the total dissipative losses 
occurring in these simulations. 

To quantify the dissipation processes in the simulations, we 
use dimensionless numbers defined as: 

- The Reynolds number, Re - SL^/v, comparing the nonlinear 
advection term to the viscous dissipation. 

- The magnetic Reynolds number, Rm - SL^/t], comparing 
magnetic field advection to the Ohmic resistivity. 

- The magnetic Prandtl number, defined as the ratio of the two 
previous quantities Pm = Rm/Re = v/rj, which measures the 
relative importance of the dissipation processes. 

In the following we use S ~ 1 as the unit of time and S L z as the 
unit of velocity. One orbit corresponds to r or b = 3nS^. 

3. Simulations of zero-net-flux MHD turbulence 

3.1. Long-timescale cycle in zero-net-flux MHD turbulence 

The first z ero-ne t-flux turbulent flow was computed by 
IStone et aD dl996l) in the context of a stratified compressible 
shearing box. In this section we consider the simpler unstrat- 
ified and incompressible case. The aspect ratio is set to L v = 
L 7 /2 and L x = 2L 7 , wh ich corresponds to the box used by 
ILesur & Longaretti! d2007l) elongated twice in the vertical direc- 
tion. As we will see, a box elongated in the z direction is useful to 



exhibit the large-scale and long-lived vertical structures one may 
miss with more classical configurations. The resolution of the 
simulation presented in this section is n x xn y xn : = 128x64x128. 
This corresponds approximately to an "equivalent" resolution 
of 256 x 128 X 256 f or a second-order finite-difference method 
dFromang et al.ll2007l) . 

We present in Fig. Q] the evolution of the azimuthal compo- 
nent of the magnetic field averaged in the x and y directions 
as a function of z and f, for a simulation with Re = 3200, 
Pm = 4. This component clearly shows some long-lived struc- 
tures, mostly of large vertical wavelength. This result is at first 
sight surprising since one expects turbulent structures to be mod- 
ified on a typical dynamical timescale, i.e. S . To quantify this 
effect more precisely, we plot in Fig. [2] the time history of the 
largest vertical Fourier mode, defined by: 

B x (k z ) = 1 f f T B x (x,y,z)exp(-ik z z)dxdydz, (10) 
L x L y L z Jo Jo Jo 

The amplitude of the largest vertical mode (with k z - 2n/L z ) 
shows clearly the cyclic behaviour already seen in Fig. [T] with 
a typical period T ~ 30 - 50 S 1 . The phase analysis shows an 
approximately constant phase during each cycle, meaning that 
these structures are not moving vertically. We have checked that 
these structures are not a transient phenomenon, integrating one 
simulation up to t — 2000 S _1 without any significant modifica- 
tion. Note that the resistive timescale for the largest-scale mode 
k z = 2n/L z is 324 S ~' in this case (Rm = 12800). For com- 
parison with other simulations, we also compute the turbulent 
transport (Fig. |2] bottom), measured by the Shakura-Sunyaev- 
like coefficient 



(B x B y ) - (v x v y ) 
S 2 Q 



(ID 



The turbulent transport obtained in our case is comparable with 
previous simulations (a ~ 7 x 10~ 3 ). However, the cyclic be- 
haviour described above is not seen in the transport coefficients 
for this simulation, which partly explains why it has not been 
observed in previous works. 

3.2. Cycle analysis 

Naturally, one may wonder what mechanism generates this mag- 
netic field structure, and whether this mechanism is related to 
some turbulent transport properties. To investigate these ques- 
tions, we first reduce the Reynolds number of the simulation, 
keeping Pm constant. This allows us to have "cleaner" flows 
to work with, removing much of the small-scale variations. We 
are able to maintain MHD turbulence and cycles down to R e = 
1600, consistent with previous results dFromang et alj2007l) . We 
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Fig. 3. Volume rendering of B x for a simulation with Re = 1600 and Km = 6400. Left panel: t = 260, corresponding to a cycle 
maximum. Right panel: t = 280, corresponding to a cycle minimum. One easily observes the large-scale B x (z) on the t = 260 
snapshot whereas strong nonaxisymmetric structures destroy the large scale structures at t — 280. 
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Fig. 4. Amplitude projected on B x (ko, f) (left) and phase (right) of the terms involved in equation (fT2l ). The amplitude clearly shows 
one cycle comparable to Fig. [2] for B x . The cycle behaviour comes mainly from the shear term, whose contribution is initially 
positive (f < 270) and later becomes negative (f > 270). EMF and resistive term always make dissipative contributions. 



show in Fig. [3] some snapshots of the azimuthal magnetic field 
B x , demonstrating the large scale field structure and its destruc- 
tion, as expected from the Fourier analysis. We then isolate one 
cycle from a Re - 1600, Pm = 4 simulation and compute the 
budget of the Fourier component B x (k z - 2n/L z ) = B x (Icq). 
Applying the transformation ( fTOb to equation ©, this budget 
may be written as: 



d,B x {ko) = S By(k ) - iko&y(k ) - r/k B x (k ), 
d,By(k ) = ik S x (ko) - rjt^B y (ko), 



(12) 
(13) 



where we have defined the electromotive force (EMF) £ = 
v X B. According to these equations, B x {ko) has three sources: 
the shear of radial magnetic field lines, the radial EMF and the 
resistivity (acting as a sink) whereas only the azimuthal EMF 
and resistivity effects appear for B Y . To study the contribution of 
each term to the global behaviour of (B x , Z? v ), we project them 



on the complex vector (B x , B^(ko), defining: 

bj = \Bj(k Q ,t)\, (14) 

nlsB y (k ,t)BAko,t)\ 

mean shear = — , (15) 

\B x (k ,t)\ 

(ko,t)Bj (ko,t)\ 

emf = '-, (16) 

\Bj(k ,t)\ 

resistive term = -rjk^B j(ko, t)\, (17) 

j being the magnetic field component studied ( j - x, y) and Sj m „ 
being the Levi-Civita tensor. These quantities, with the phases 
of their associated terms in (fT2l -(fT3l). are plotted in Fig. H]for 
budget ( fT2l ) and in Fig. [5] for budget ( fT3l . 

We find in Fig. [4] the long-timescale cycle we have already 
described. The shear term is positive in the beginning of the cy- 
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Fig. 5. Amplitude projected on B y (ko, f) (left) and phase (right) of the terms involved in equation ( TT3b - The amplitude of B y is 
approximately 20 times smaller than B x . The cycle observed for B Y is related to the oscillation of the EMF £ v which is reversed at 
t 243. 



cle, but becomes negative for t > 270, whereas the EMF has 
a systematic resistive effect, with a clear phase correlation be- 
tween the EMF and the resistive term. According to these re- 
sults, the long-timescale cycle in the azimuthal field comes from 
the behaviour of the radial field, whereas the radial EMF may 
be seen, in first approximation, as a turbulent resistivity. As ex- 
pected from this conclusion, B v is also oscillating on a long 
timescale (Fig. |5j with a smaller amplitude compared to B x . 
Moreover, this cycle comes clearly from an oscillation of the 
azimuthal EMF <5 A , which is reversed at t ~ 243. Note, however, 
that the EMF is very small (~ 5 x 10~ 4 ) compared to quantities 
in the budget for the azimuthal field, which may explain the long 
period of these oscillations. 

3.3. Modal analysis of non axisymmetric structures 

The EMFs described previously are obviously nonlinear terms, 
involving a coupling between v and B. An interesting question 
is therefore which modes contribute most to the EMFs observed 
in Figs [4^5] In particular, one may wonder whether nonaxisym- 
metric structures play a role and if so, which structures are dom- 
inant. Following this idea, we compute the following quantities: 

£(»,,£<)) = tV f T *2K[v(n x ,y,z) X B(n x ,y,zT] 

L y L z Jo Jo 

x exp(-/fc z) dy dz, (18) 

with the fields Tp(n x , y, z) defined by the Fourier coefficient in the 
x direction: 

V(n x ,y,z) = J (p(x,y,z)exp^-2im x j-^dx. (19) 

In this notation, £>{n x , ko) corresponds to the contribution of the 
nonaxisymmetric mode n x to the total EMF G(ko). We then put 
G(n x , ko) in equation ( [ToT l and plot the contribution of the first 
few nonaxisymmetric modes in Fig. [6] From this figure, it is 
clear that most of the EMF comes from the largest azimuthal 
wavelengths, n x < 3. However, the detailed generation of S x and 
S y is somewhat different: although the contribution of axisym- 
metric modes (n x = 0) to S x is small, these modes are clearly 
important for £ y . Therefore, the azimuthal EMF S x comes es- 
sentially from the largest nonaxisymmetric mode whereas the 
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Fig. 7. Magnetic helicity history during one cycle. The plain 
curve is the helicity due to the large-scale field k z = 2n/L z , 
the dashed line corresponds to the helicity associated with the 
largest nonaxisymmetric modes, and the total helicity is plotted 
as a dash-dotted line. 



radial EMF comes from a more general effect of both nonax- 
isymmetric and axisymmetric structures. These results show that 
the magnetic cycles observed in the simulations are primarily 
due to three-dimensional structur es, as one wo uld expect from 
Cowling's antidynamo theorem (Cowling 198~I1). 

It has been pointed out by Blac kmanl (120071) that monitoring 
the helicity evolution may be key to understanding large-scale 
field production. Moreover, since £• B is non-zero for the largest 
vertical modes, some large-scale magnetic helicity should be 
produced. To explore the helicity behaviour in our case, we com- 
pare the box-averaged magnetic helicity due to the vertical mode 
k = (2n/L z )e z with the magnetic helicity associated with the 
largest nonaxisymmetric modes n$ = 1,2 in Fig. [Vj During one 
cycle, we note that the helicity associated with the vertical modes 
changes sign when B v is reversed (f 270). Moreover, this rever- 
sal seems to be associated with an exchange of magnetic helic- 
ity between the nonaxisymmetric and the axisymmetric vertical 
modes. Note however that the total magnetic helicity also varies 
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significantly during one cycle. Since the box averaged total mag- 
netic helicity is conserved in ideal MHD, any variation of the to- 
tal helicity must be due to a resistive (and therefore small-scale) 
process. In that way, we cannot make a very clear distinction be- 
tween the helicity exchanged with the nonaxisymmetric modes, 
and the helicity due to small scale production without any further 
analysis. Moreover, we emphasise that the averaged helicity ap- 
pears to be about 10~ 2 times smaller than L z (b 2 ), which is prob- 
ably a consequence of the numerous symmetries in unstratified 
shearing boxes. Therefore, the helicity produced in our simula- 
tions might be more a consequence of the underlying nonlinear 
dynamo process rather than its cause. 

In this section, we have shown that zero-net-flux MRI simu- 
lations exhibit long-timescale oscillations in the azimuthal mag- 
netic field B x with a large vertical wavelength. Studying the bud- 
get of B x , it appears that these cycles are primarily an ampli- 
fication of oscillations in B v through the shear term, whereas 
the EMF £, acts as a turbulent resistivity. A similar analysis of 
the By cycles shows that they are generated by a long-period 
but small-amplitude oscillation of £ v . Finally, we used a modal 
analysis to demonstrate that £ x and £,, are mostly generated by 
large-scale nonaxisymmetric waves, making these cycles a fully 
three-dimensional problem. 

Therefore, to have a better understanding of this phenomeno- 
logical picture, one needs to study these nonaxisymmetric struc- 
tures and the EMFs associated with them, which is the topic of 
the next section. 



4. Linear analysis of nonaxisymmetric waves 

4.1. Model and equations 

To understand the nonaxisymmetric origin of the EMF described 
previously, we consider a linear model including shearing waves 
in the presence of a background azimuthal magnetic field with 
a vertical structure, B^(z). This vertical structure is required in 
order to compare the EMF generated by the perturbation with 
the background magnetic field. In this sense, this calculation dif- 
fers from the initia l MRI studies with a uniform azimuthal field 
dBalbus & HawlevlfT992h . 

Because of the mean shear, we cannot use a classical Fourier 
decom position for nonaxisymmetric w ave s. Therefore, fol- 
lowing iGoldreich & Lynden-Belll d 19651) and iBalbus & Hawleyl 



(119921) . we define a sheared frame comoving with the laminar 
flow: 

x' = x - Syt y' =y z' = z t' = t. (20) 

In this Lagrangian frame, we look for plane-wave solutions <f> = 
cp(z' , t') exp[i(k' x x' + k' y y')] which may be transformed into the 
unsheared frame as: 

(p(x,y,z,t) = <p(z,f)exp[i(k x x + k y (t)y)], (21) 

with k x = k' x and k y (t) = kf y — Sk x t. We then consider per- 
turbations from the background magnetic structure in the form 
of these shearing waves. This vertical structure is supposed to 
mimic the large-scale field observed in numerical simulations, 
i.e. being a function of z only. Moreover, since the cycles have a 
long timescale compared to the shear, we also assume to a first 
approximation that this large-scale field is constant for a shear- 
ing wave. We therefore write the general linearized solution as: 

v = v(z, t) exp[i(k x x + ky(t)y)], (22) 
B = b{z, t) exp[i(k x x + k y (t)y)] + B° x (z)e x , (23) 

where \v\ and |b| are supposed infinitely small. Using these solu- 
tions in the evolution equations (H)-©, one eventually finds: 

d,v = -(ik + e z d z )fl + (2Q - S)v y e x - 2D.v x e y 

+ik x B° x (z)b + b z d z B° x (z)e x + v(d 2 - k 2 )v, (24) 
d,b = Sb y e x + ik x B x (z)v - v z d z B x (z)e x + T](d 2 z - k 2 )b, (25) 

where the generalised pressure is constrained by the incompress- 
ibility condition {7]i and k 2 = k 2 + k 2 (t). 

4.2. Numerical solution 

In the following, we will assume a large-scale field similar to 
the oscillating Fourier mode studied in the previous section. We 
therefore set: 

B° x (z) = B cos(A:oz). (26) 

The resistive diffusion of this non-uniform field may be ne- 
glected on the timescales of interest here. To solve the equa- 
tions (l24l and ( |25] l, we consider a shearing wave (k x , k y {t)), L z - 
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Fig. 8. Evolution of T x (upper panels) and F v (lower panels) for a set of 40 shearing waves with random initial conditions, Re = 1600 
and Rm = 6400. From left to right, the large-scale field is increased with Bq = 0.02; 0.04; 0.08; 0.2. (F A ) is always negative whereas 
(T y ) is positive for Bq < 0.08 and is reversed for larger field strengths. This is consistent with the cycles observed in the previous 
section (see text). 
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Fig. 9. Evolution of T x (upper panels) and F, (lower panels) for a set of 40 shearing waves with random initial conditions, Re = 6400 
and Rm = 25600. From left to right, the large-scale field is increased with Bq = 0.04; 0.08; 0.2; 0.3. As for Fig. [8] (T x ) is always 
negative whereas (Fy) is reversed in this case for Bq > 0.2. 



periodic in the vertical directiorQ. Without any loss of general- 
ity, we also define a new time variable fg W so that feytew) - 0. 
To check whether the large-scale field has a systematic impact 
on the structure of shearing waves, we compute numerically the 
evolution of several shearing waves using random initial condi- 
tions. The shearing waves are initialised at fsw = -10 with ran- 
dom phases and randomly oriented velocity and magnetic fields. 
The initial velocity and magnetic fields are verified to satisfy (Q 
and ([8]). We then evolve these shearing waves using a Fourier de- 
composition with n, = 64 modes in the vertical direction up to 
t = 30. Finally, we compute the axisymmetric EMF due to each 
shearing wave, £ sw (z) = 2%{v x b*), and check the correlation 
of its curl with the imposed large-scale field, i.e. 

J Q L: B° x (z)d z S™(z)dz 

Fy — Sjzm D • (27) 



1 The periodicity in the vertical direction is required since the source 
term B°(z) is periodic with period L : . 



We plot r x and F v for a set of 40 shearing waves with Re = 1600, 
Rm = 6400, and the aspect ratio used in the last section, in Fig. [8] 
We have considered only the largest nonaxisymmetric waves 
k x - 2n/L x as they appear to be the dominant ones in the non- 
linear simulations (see section [331 ). Although the F coefficients 
are assumed to be infinitely small compared to Bq, they provide 
an understanding of how a finite-amplitude shearing wave may 
modify the large-scale field B° K (z). To clarify this idea, let us now 
assume that B Q x (z) can vary slowly in time. Then, using (fT2l and 
( fT3l , one find that T x > corresponds to a increase of the large- 
scale field B° x whereas T x < is equivalent to a resistive effect. 
Following the same argument, F v > corresponds to a creation 
of a large-scale B v (z) with the same sign as B x (z) which can, in 
combination with the shear term, amplify B x (z). On the other 
hand, F v < leads to a B, (z) with the opposite sign to B° K (z) and 
possibly, in combination with the shear, to a destruction of B° K . 
Note however that we have not discussed here the amplitudes 
of the T coefficients but just their signs. The main reason is that 
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Fig. 10. Evolution of (b x b y ) (upper panels) and (v x v y ) (lower panels) for a set of 40 shearing waves with random initial conditions, 
Re = 6400 and Rm = 25600. From left to right, the large-scale field is increased with Bq = 0.04; 0.08; 0.2; 0.3. As expected from 
traditional MRI results, (b x b y ) > and (v^y,) < 0: the shearing waves transport angular momentum outwards. 



their final amplitude will depend quadratically on the initial ex- 
citation of the shearing waves. In a real system, this excitation 
is a highly nonlinear process which depends on the small-scale 
properties of the turbulence^ Therefore, our random excitation 
is not a precise enough turbulence model and we cannot rely on 
the amplitude found in this linear analysis. Furthermore, the evo- 
lution of the shearing waves may differ from the predictions of 
this analysis if they reach nonlinear amplitudes. 

As one can see in Fig. [8] V x has a negative sign on average, 
independent of Bo- As mentioned previously, this can be inter- 
preted as a resistive effect on B x and is similar to the systematic 
resistive effect observed in numerical simulations (see Fig. |4|. 
This effect is easily explained as resulting from the mixing of the 
non-uniform azimuthal field by vertical motions. The behaviour 
of r v is a little more complicated since it reverses for Bq ~ 0.08. 
For small Bo, F,, > which implies a possible growth of B y (z) 
and therefore an increasing positive shear term in (fT2l . On the 
contrary, for Bo > 0.08, one expects the shear term to decay and 
eventually be reversed. This property is indeed observed in nu- 
merical simulations. In particular, the shear term in Fig.|4]starts 
to decay for Bo > 0.08 which corresponds to the the predicted 
result from the linear analysis. 

To check the effect of dissipation processes on this picture, 
we plot in Fig.|9]the evolution of the F coefficients for Reynolds 
and magnetic Reynolds numbers 4 times larger than in the pre- 
vious case. One finds essentially the same properties: a system- 
atic resistive effect for T x and a reversal for F v . Note however 
that the reversal appears for a larger Bo in the less diffusive case 
(Bq ~ 0.2). We have also checked that these properties persist 
when one changes the magnetic Prandtl number at sufficiently 
large Re and Rm. These remarks suggest that the results found in 
this linear analysis are not related to a finite resistivity or vis- 
cosity property and may therefore appear for arbitrarily large 
Reynolds numbers. 



2 Note that these properties may depend in turn on the large-scale 
field and the amplitudes of the shearing waves. 



4.3. Phenomenological properties 

The waves described in this section are clearly inhomoge- 
neous in the vertical direction. However, we can understand 
them as a version of the magnetorot ational instability in the 
presence of a varying azimut hal field (Bal bus & HawleyH l992; 
Terq uem & PapaloizouH l996). As one would expect, the trans- 
port coefficients (b x b y ) and (v x v y ) associated with the linear 
waves are nonzero, and lead to an outward angular momentum 
trans port (Fig. [10|, as in the classical vertical field calculation 
(see iPessah et alfl2008l) . This demonstrates that these shearing 
waves actually extract energy from the mean flow as a classi- 
cal MRI mode. Because of the resistive properties of S y , these 
waves also extract some energy from the large-scale azimuthal 
field, which is generated through shearing of the radial field. 
Therefore, the energy of the waves comes primarily from the 
mean shear, as expected. Moreover, it is known that the az- 
imuthal MRI has an optimum growth rate for a given azimuthal 
wavelength k x and field strength B x . In the high-A:, limit with a 
uniform B - Bpe^, this maximum growth rate is reached when 
k x Bo = V5/125 (lOgilvie & Pringlelll996l) . Using our parame- 
ters, this leads to Bo ~ 0.25 L z , which is surprisingly close to the 
magnetic field amplitude for which F x is reversed at large Re. 
Although this argument is by no mean comprehensive, it draws 
attention to the close relation between the transient amplifica- 
tio n observed in our analysis and the (ideal) MRI waves studied 
by iBalbus & Hawlevl (Q~992). Naturally, these conjectures need 
to be checked more carefully using a proper analytical analysis, 
which will be the subject of another paper. 

5. A toy model 

In this section, we provide a toy model reproducing the basic lin- 
ear properties exhibited in the previous section. This toy model 
does not pretend to be an accurate set of closure relations for 
equations (fT2l-(fT3l but it includes the main physical ingredients 
required to reproduce qualitatively the cycle behaviour described 
in section [3] We therefore rewrite equations (fT2l-([T3l> as: 

d t B x (k , t) = S B y (ko, t) - BB x (k , t-t r ), (28) 

a on a nn t . ,B r - \B x (ko,t - t r )\ 

d,By(ko,t) = yB x (koJ-t r ) , (29) 
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Fig. 11. Numerical solution of the toy model (equations |28T|29T > 
with y = 0.007,/? = 0.010, t r = 2 5"' and B r = 0.08, exhibit- 
ing sustained magnetic field oscillations. The time between two 
B x reversals is about ~ 50 S^ 1 as in the full simulations. This 
toy model reproduces qualitatively the results of our numerical 
simulations. 



quantities. Interestingly, the cycles obtained using the model are 
self-sustained despite the presence of a resistive effect, showing 
that the basic properties discussed previously are sufficient to in- 
ject magnetic energy into the system. One should note however 
that the typical EMF term involved in the evolution of B x is sig- 
nificantly smaller in the model than in the simulation. Moreover, 
we observe a perfect field reversal of B x which is clearly not 
present in the simulations. This may be explained by the pres- 
ence of "noise" due to small turbulent eddies which are continu- 
ously exciting the B x (ko) mode in the simulations. 

As we have seen, our toy model is able to reproduce the ba- 
sic properties observed in complicated nonlinear simulations, 
just considering the linear results from the previous section. 
However, this model is far from perfect. For example, lineariz- 
ing the model around B = leads to the wrong conclusion that 
the trivial solution B — is unstable. This peculiar property 
comes from the assumption of a constant turbulent background. 
Therefore, if we assume that something (e.g. turbulence) is con- 
stantly exciting shearing waves, then the solution B — is unsta- 
ble. However, in a real system, the turbulent motions are partly 
generated by these shearing waves when they reach a nonlin- 
ear regime. If the amplification is not large enough (i.e. around 
B = 0), the shearing waves cannot sustain the turbulence, and 
therefore the flow cannot excite new shearing waves: the whole 
system relaxes to B = 0. In conclusion, despite the apparent 
simplicity of our toy model, the laminar flow is not subject to 
a linear instability and the underlying mechanism sustaining the 
turbulence in this system is clearly a nonlinear effect, potentially 
involving the transient amplification described earlier. 



where y,B and B r are constants. In this set of equations, we have 
neglected the physical dissipation processes, as they are very 
small in the budget of the simulations (see Figs |4]-[5j- We have 
assumed a resistive term for 8> y , as expected from the linear anal- 
ysis. However, to take into account the fact that the EMF is due 
to a progressively amplified shearing wave, we assume the EMF 
is slightly retarded with respect to the large-scale magnetic field, 
with a delay time t r . Since the typical shearing wave growth rate 
is of the order of a shear time, we expect in first approximation 
t r ~ S . The model used for <5 A reproduces the reversal de- 
scribed in the linear analysis, at \B X \ = B r . It also includes the 
delay used in the £ v model for the same reasons. This term will 
be referred to in the following as the y effect term, as it is not 
formally equivalent to the classical a effect used in dynamo the- 
ory, but more an effect of anisotropic turbulent resistivity. Note 
that the same kind of model i nvolving an anisotropic tu r bulen t 
resistivity has been studied by Rogachevski i & Kleeorinl f2003) 
in the context of the shear dynamo. For this model to be con- 
sistent with the previous analysis, we have to assume that the 
underlying three-dimensional flow is turbulent in some way, so 
that small-amplitude shearing waves are continuously excited 
for fsw < 0, and then amplified linearly. Therefore, this model 
assumes that the flow is already subject to a three-dimensional 
turbulence and describes the variations of the large-scale field. 

We numerically solve the toy model using y = 7 X 10~ 3 , 
B = 10~ 2 , t r = 2S~ 1 and B r = 8 X 10~ 2 . Here y and B are chosen 
so that our model mimics the main behaviour of the numerical 
simulations, whereas B, is chosen according to the results of the 
linear analysis. The resulting evolution of B x and B y is plotted 
in Fig.[TT] where the "mean shear" curve corresponds to the first 
term on the right-hand side of equation ((28), and the EMFs are 
the B and y terms. When comparing with the fully nonlinear cy- 
cle (Figs |4]-[5]l, we find essentially the same time history for all 



6. Discussion 

6.1. Summary 

In this paper, we have investigated the behaviour of the large- 
scale magnetic field in zero-net-flux simulations of the magne- 
torotational instability. We have first shown that the large-scale 
azimuthal field B x (z) is subject to a long-timescale oscillation 
when the flow is turbulent. Studying the induction equation, we 
have seen that these cycles are primarily due to an oscillation 
of a large-scale radial magnetic field B y and an azimuthal EMF 
& x , whereas the radial EMF <5 V acts like a turbulent resistivity 
on B x . Studying the nonaxisymmetric modes of the system, we 
found that the axisymmetric EMF is due mostly to the largest 
nonaxisymmetric structures, showing that the cycles are essen- 
tially a large-scale process. To understand the generation of the 
EMF, we have investigated the behaviour of nonaxisymmetric 
linear waves in the presence of an inhomogeneous and constant 
in time azimuthal magnetic field B x (z). This linear analysis ex- 
plains most of the properties of the EMF observed in the nonlin- 
ear simulations. In particular, we have shown that the large-scale 
axisymmetric S x can be reversed for large B x , explaining the 
cycles in the simulations. To summarise these results, we have 
considered a simplified closure model, encapsulating the linear 
properties. This model is able to reproduce the main behaviours 
of the cycles despite its simplicity, showing that the physics in- 
volved in the cycles is well described by our linear analysis. 

To summarise our findings, we sketch the mechanism re- 
sponsible for the cycles in Fig. [12] In this sketch, we have 
assumed that the nonaxisymmetric excitation (seed) is related 
to the nonlinear coupling of amplified shearing waves (dashed 
line). This specific process has not been studied in the present 
work and is just given here as a way to "close the loop". 
Note however that a more complicated mechanism may be in- 
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Fig. 12. Scheme of the mechanism responsible for the large- 
scale magnetic field cycles. The dashed lines represent the non- 
linear interactions which are not discussed in the present work. 
These interactions are nonetheless required to obtain a self- 
sustained mechanism (see text). 



volved in this back-reaction without changing the global pic- 
ture. Moreover, we have described the contribution of & x (z) as a 
generic "dynamo" effect, but one should remember that this term 
may be either correlated or anticorrelated with B x , as a function 
of the strength of B x . Finally, we have omitted viscous and re- 
sistive effects for simplicity, as they have a small impact on the 
mechanism itself. 

6.2. Comparison with previous works 

We would like to stress that the model presented here does not 
constitute a full description of a sustaining mechanism for MHD 
turbulence in discs. Indeed, we have assumed in the linear anal- 
ysis and in the closure model that the flow is able to gener- 
ate continuously small-amplitude shearing waves. As mentioned 
previously, this process involves a nonlinear coupling of non- 
axisymmetric waves, which is beyond the scope of this work. 
However, a fully consistent self-sustaining pro cess would re- 
quire a description of this effect. According to iFromang et al.l 
(120071) . MHD turbulence appears in zero-net-flux shearing-box 
simulations only when Pm > 1 . However, most of our analysis 
does not depend on the magnetic Prandtl number, and the dy- 
namo process described here may work for arbitrarily low Pm. 
This suggests that the shearing-wave excitation mechanism de- 
pends on the Prandtl number, and may be too weak in the Pm < 1 
case. This confirms that the excitation should be related to the 
small turbulent scales, which are still poorly understood. 

Interestingly, our results indicate that the large-scale field 
generation mechanism is not destroyed for large Reynolds num- 
bers, provided that the background flow is turbulent. Therefore, 
the same kind of mechanism might be at work in real astrophys- 
ical discs, generating a large-scale field in a sufficiently ionised 
and turbulent plasma. One should note however that the config- 
uration used in our simulations is not comparable to the real ge- 
ometry of an accretion disc. In particular, we do not know how 
the aspect ratio will modify the mechanism of the cycles, nor 
how the vertical stratification may enter t his picture. However, 
as discovered by Brandenburg et al. ( 1995) in the case of a strat- 
ified flow with nonperiodic vertical boundary conditions, one ob- 
serves azimuthal field cycles with a long timescale. These cycles 
have a significantly longer period (T ~ 200 S _1 ) compared to 
ours and involve a modification of the net azimuthal magnetic 
flux. M oreover, it has been shown by iBrandenburg & Donnerl 
d 19971) that in this case, the cyclic behaviour may be explained 
by a classical a effect proportional to the vorticit^Q Because 

3 Note that this a effect might be related to the presence of vertical 
stratification and the breaking of the reflectional symmetry 
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of these fundamental differences, we cannot conclude that the 
mechanism responsible for all these cycles is the same without 
further investigation. Nevertheless, it seems that a more system- 
atic investigation of this kind of long-timescale behaviour is re- 
quired in order to better understand the way MHD turbulence 

sustains magnetic fields in various con figurations. 

As mentioned in the introduction, Pess ah etldl ([2007) sug- 
gested that turbulent transport should vanish for large enough 
resolution and Reynolds number. Therefore, the turbulent trans- 
port associated with this mechanism is another important issue. 
As one can check easily in the numerical simulation and in the 
toy model, the large-scale B x and B y are roughly phase-shifted 
by 7r/2, leading to a very small Maxwell stress ((B x B y ) ~ 10~ 4 ). 
Therefore, in our picture, most of the transport must come from 
the shearing waves, which have been shown to have nonzero 
Maxwell and Reynolds stresses. Unfortunately, we cannot put 
a precise value on the transport due to the shearing waves since 
it is controlled by their initial excitation. Therefore, this mecha- 
nism cannot give a precise answer on the expected transport in 
a real accretion disc. From a more phenomenological point of 
view, we expect the large-scale magnetic field strength produced 
by our mechanism not to depend on the Reynolds numbers if 
they are sufficiently large (see section l4~2l >. Since this large-scale 
field couples with all the other modes available, we expect the 
large-scale nonaxisymmetric structures to have roughly the same 
amplitude as the large-scale field, in the limit of a fully devel- 
oped turbulence. If this picture is correct, we would then expect 
a minimum transport of the order of a few times 10~ 3 , indepen- 
dently of the Reynolds numbers. Note however that this conclu- 
sion relies on the assumption that the large-scale field does not 
depend on the dissipation coefficients. We have shown that this 
assumption is plausible, although a precise analytical descrip- 
tion of our linear analysis would be required to ascertain this 
statement. 

The present description of an accretion disc dynamo is 
related to the stea dy self-sustaining solutions obtained by 
iRincon et al.l d2007l) by numerical continuation methods. Both 
mechanisms involve a large-scale azimuthal magnetic field with 
zero net flux, which is generated through the shearing of a radial 
field. Both also involve a nonaxisymmetric magnetorotational- 
type instability that regenerates the radial field through its non- 
linear feedback. However, the present mechanism is intrinsically 
local and unsteady, and may work for arbitrarily large Reynolds 
numbers whereas t he steady, highly symmetrical solutions of 
IRincon et al.l d2007) depend on the existence of walls and are ap- 
parently restricted to low Reynolds numbers. Nevertheless, one 
may also think that these steady solutions correspond to some 
fixed points in phase space, around which the shearing box solu- 
tions oscillate, creating the dynamo cycles we have described. If 
this picture is correct, it would be an elegant way to unify these 
different approaches, and possibly to find other mechanisms of 
the same kind. 

This w ork might also be c ompared to the shear dynamo sim- 
ulations of lYousef et al.l(l2008l) . However, we emphasise that our 
analysis involves a nonlinear dynamo, whereas the shear dy- 
namo described by lYousef et al.l {2008) is a kinematic (linear) 
dynamo, in which turbulence is forced. Moreover, the shear dy- 
namo occurs in a non-rotating system in which the MRI does not 
appear. Applying our linear analysis to such a system does not 
produce a y effect that could explain the shear dynamo in either 
linear or nonlinear regimes. 

Finally, we would like to stress that the problem of a non- 
linear dynamo, in volving several lin ear instabilities, has also 
been described by ICline et ail d2003l) in the context of magnet- 
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ically buoyant flows. Although the underlying instabilities are 
somewhat different, they also found a cyclic behaviour which, 
as in the present work, cannot be fully described by a classical 
a effect. This may indicate that all these instabilities share some 
common properties yet to be discovered, allowing for the devel- 
opment of a nonlinear dynamo. 

6.3. Future work 

As discussed previously, the main issue raised by our find- 
ings is how the turbulence is able to excite shearing waves. 
Interestingly, the same kind of problem ari ses in the case of 
the sh earing box with a mean azimuthal field (B albus & Hawlevl 
1 1 992b . Therefore, a simpler way to study this effect is to in- 
vestigate the way turbulence is sustained in the presence of a 
mean azimuthal field. In particular, an interesting test would 
be to check that in this case, the prese nce of turbulence als o 
depends on Pm, in a similar way as in Froma ng et al.l (|2007). 
This would be a good indication that the large-scale dynamo 
process itself does not depend on Pm at large Re, even though 
the underlying turbulent regeneration mechanism process does. 
Moreover, this idea would provide a new way to understand the 
fm-dependence obs erved in MRI simulations with a nonzero 
mean vertical field (Lesur & Longaretti 2007), and potentially 
to describe the asymptotic transport in the limit of high and low 
Pm. 

On the other hand, the dynamo process itself requires fur- 
ther investigation. The linear analysis provided here has been 
computed using essentially numerical tools, as the analytical de- 
scription of this problem is somewhat technical. This analytical 
description is nevertheless required, as it will give a precise and 
formal constraint on the y effect and its dependence on the az- 
imuthal field strength. It may also be a way to have a more phys- 
ical understanding of the underlying mechanism responsible for 
the resistive and y effects, leading to a potential generalisation 
of this mechanism to a wider class of flow. Finally, we would be 
able to check, at least linearly, if this dynamo process depends on 
the non-ideal effects, and potentially to confirm our conjecture. 
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